home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / EXAMPLE.CPP < prev    next >
C/C++ Source or Header  |  1995-01-11  |  10KB  |  281 lines
  1. //$$ example.cpp                             Example of use of matrix package
  2.  
  3. #define WANT_STREAM                  // include.h will get stream fns
  4. #define WANT_MATH                    // include.h will get math fns
  5.                                      // newmatap.h will get include.h
  6.  
  7. #include "newmatap.h"                // need matrix applications
  8.  
  9. #include "newmatio.h"                // need matrix output routines
  10.  
  11.  
  12. // demonstration of matrix package on linear regression problem
  13.  
  14.  
  15. void test1(Real* y, Real* x1, Real* x2, int nobs, int npred)
  16. {
  17.    cout << "\n\nTest 1 - traditional\n";
  18.  
  19.    // traditional sum of squares and products method of calculation
  20.    // with subtraction of means
  21.  
  22.    // make matrix of predictor values
  23.    Matrix X(nobs,npred);
  24.  
  25.    // load x1 and x2 into X
  26.    //    [use << rather than = with submatrices and/or loading arrays]
  27.    X.Column(1) << x1;  X.Column(2) << x2;
  28.  
  29.    // vector of Y values
  30.    ColumnVector Y(nobs); Y << y;
  31.  
  32.    // make vector of 1s
  33.    ColumnVector Ones(nobs); Ones = 1.0;
  34.  
  35.    // calculate means (averages) of x1 and x2 [ .t() takes transpose]
  36.    RowVector M = Ones.t() * X / nobs;
  37.  
  38.    // and subtract means from x1 and x1
  39.    Matrix XC(nobs,npred);
  40.    XC = X - Ones * M;
  41.  
  42.    // do the same to Y [use Sum to get sum of elements]
  43.    ColumnVector YC(nobs);
  44.    Real m = Sum(Y) / nobs;  YC = Y - Ones * m;
  45.  
  46.    // form sum of squares and product matrix
  47.    //    [use << rather than = for copying Matrix into SymmetricMatrix]
  48.    SymmetricMatrix SSQ; SSQ << XC.t() * XC;
  49.  
  50.    // calculate estimate
  51.    //    [bracket last two terms to force this multiplication first]
  52.    //    [ .i() means inverse, but inverse is not explicity calculated]
  53.    ColumnVector A = SSQ.i() * (XC.t() * YC);
  54.  
  55.    // calculate estimate of constant term
  56.    //    [AsScalar converts 1x1 matrix to Real]
  57.    Real a = m - (M * A).AsScalar();
  58.  
  59.    // Get variances of estimates from diagonal elements of invoice of SSQ
  60.    //    [ we are taking inverse of SSQ - we need it for finding D ]
  61.    Matrix ISSQ = SSQ.i(); DiagonalMatrix D; D << ISSQ;
  62.    ColumnVector V = D.AsColumn();
  63.    Real v = 1.0/nobs + (M * ISSQ * M.t()).AsScalar();
  64.                         // for calc variance of const
  65.  
  66.    // Calculate fitted values and residuals
  67.    int npred1 = npred+1;
  68.    ColumnVector Fitted = X * A + a;
  69.    ColumnVector Residual = Y - Fitted;
  70.    Real ResVar = Residual.SumSquare() / (nobs-npred1);
  71.  
  72.    // Get diagonals of Hat matrix (an expensive way of doing this)
  73.    Matrix X1(nobs,npred1); X1.Column(1)<<Ones; X1.Columns(2,npred1)<<X;
  74.    DiagonalMatrix Hat;  Hat << X1 * (X1.t() * X1).i() * X1.t();
  75.  
  76.    // print out answers
  77.    cout << "\nEstimates and their standard errors\n\n";
  78.    cout.setf(ios::fixed, ios::floatfield);
  79.    cout << setw(11) << setprecision(5)  << a << " ";
  80.    cout << setw(11) << setprecision(5)  << sqrt(v*ResVar) << endl;
  81.    // make vector of standard errors
  82.    ColumnVector SE(npred);
  83.    for (int i=1; i<=npred; i++) SE(i) = sqrt(V(i)*ResVar);
  84.    // use concatenation function to form matrix and use matrix print
  85.    // to get two columns
  86.    cout << setw(11) << setprecision(5) << (A | SE) << endl;
  87.    cout << "\nObservations, fitted value, residual value, hat value\n";
  88.    cout << setw(9) << setprecision(3) <<
  89.      (X | Y | Fitted | Residual | Hat.AsColumn());
  90.    cout << "\n\n";
  91. }
  92.  
  93. void test2(Real* y, Real* x1, Real* x2, int nobs, int npred)
  94. {
  95.    cout << "\n\nTest 2 - Cholesky\n";
  96.  
  97.    // traditional sum of squares and products method of calculation
  98.    // with subtraction of means - using Cholesky decomposition
  99.  
  100.    Matrix X(nobs,npred);
  101.    X.Column(1) << x1;  X.Column(2) << x2;
  102.    ColumnVector Y(nobs); Y << y;
  103.    ColumnVector Ones(nobs); Ones = 1.0;
  104.    RowVector M = Ones.t() * X / nobs;
  105.    Matrix XC(nobs,npred);
  106.    XC = X - Ones * M;
  107.    ColumnVector YC(nobs);
  108.    Real m = Sum(Y) / nobs;  YC = Y - Ones * m;
  109.    SymmetricMatrix SSQ; SSQ << XC.t() * XC;
  110.  
  111.    // Cholesky decomposition of SSQ
  112.    LowerTriangularMatrix L = Cholesky(SSQ);
  113.  
  114.    // calculate estimate
  115.    ColumnVector A = L.t().i() * (L.i() * (XC.t() * YC));
  116.  
  117.    // calculate estimate of constant term
  118.    Real a = m - (M * A).AsScalar();
  119.  
  120.    // Get variances of estimates from diagonal elements of invoice of SSQ
  121.    DiagonalMatrix D; D << L.t().i() * L.i();
  122.    ColumnVector V = D.AsColumn();
  123.    Real v = 1.0/nobs + (L.i() * M.t()).SumSquare();
  124.  
  125.    // Calculate fitted values and residuals
  126.    int npred1 = npred+1;
  127.    ColumnVector Fitted = X * A + a;
  128.    ColumnVector Residual = Y - Fitted;
  129.    Real ResVar = Residual.SumSquare() / (nobs-npred1);
  130.  
  131.    // Get diagonals of Hat matrix (an expensive way of doing this)
  132.    Matrix X1(nobs,npred1); X1.Column(1)<<Ones; X1.Columns(2,npred1)<<X;
  133.    DiagonalMatrix Hat;  Hat << X1 * (X1.t() * X1).i() * X1.t();
  134.  
  135.    // print out answers
  136.    cout << "\nEstimates and their standard errors\n\n";
  137.    cout.setf(ios::fixed, ios::floatfield);
  138.    cout << setw(11) << setprecision(5)  << a << " ";
  139.    cout << setw(11) << setprecision(5)  << sqrt(v*ResVar) << endl;
  140.    ColumnVector SE(npred);
  141.    for (int i=1; i<=npred; i++) SE(i) = sqrt(V(i)*ResVar);
  142.    cout << setw(11) << setprecision(5) << (A | SE) << endl;
  143.    cout << "\nObservations, fitted value, residual value, hat value\n";
  144.    cout << setw(9) << setprecision(3) <<
  145.       (X | Y | Fitted | Residual | Hat.AsColumn());
  146.    cout << "\n\n";
  147. }
  148.  
  149. void test3(Real* y, Real* x1, Real* x2, int nobs, int npred)
  150. {
  151.    cout << "\n\nTest 3 - QR triangularisation\n";
  152.  
  153.    // QR triangularisation method
  154.  
  155.    // load data - 1s into col 1 of matrix
  156.    int npred1 = npred+1;
  157.    Matrix X(nobs,npred1); ColumnVector Y(nobs);
  158.    X.Column(1) = 1.0;  X.Column(2) << x1;  X.Column(3) << x2;  Y << y;
  159.  
  160.    // do Householder triangularisation
  161.    // no need to deal with constant term separately
  162.    Matrix X1 = X;                 // Want copy of matrix
  163.    ColumnVector Y1 = Y;
  164.    UpperTriangularMatrix U; ColumnVector M;
  165.    QRZ(X1, U); QRZ(X1, Y1, M);    // Y1 now contains resids
  166.    ColumnVector A = U.i() * M;
  167.    ColumnVector Fitted = X * A;
  168.    Real ResVar = Y1.SumSquare() / (nobs-npred1);
  169.  
  170.    // get variances of estimates
  171.    U = U.i(); DiagonalMatrix D; D << U * U.t();
  172.  
  173.    // Get diagonals of Hat matrix
  174.    DiagonalMatrix Hat;  Hat << X1 * X1.t();
  175.  
  176.    // print out answers
  177.    cout << "\nEstimates and their standard errors\n\n";
  178.    ColumnVector SE(npred1);
  179.    for (int i=1; i<=npred1; i++) SE(i) = sqrt(D(i)*ResVar);
  180.    cout << setw(11) << setprecision(5) << (A | SE) << endl;
  181.    cout << "\nObservations, fitted value, residual value, hat value\n";
  182.    cout << setw(9) << setprecision(3) << 
  183.       (X.Columns(2,3) | Y | Fitted | Y1 | Hat.AsColumn());
  184.    cout << "\n\n";
  185. }
  186.  
  187. void test4(Real* y, Real* x1, Real* x2, int nobs, int npred)
  188. {
  189.    cout << "\n\nTest 4 - singular value\n";
  190.  
  191.    // Singular value decomposition method
  192.  
  193.    // load data - 1s into col 1 of matrix
  194.    int npred1 = npred+1;
  195.    Matrix X(nobs,npred1); ColumnVector Y(nobs);
  196.    X.Column(1) = 1.0;  X.Column(2) << x1;  X.Column(3) << x2;  Y << y;
  197.  
  198.    // do SVD
  199.    Matrix U, V; DiagonalMatrix D;
  200.    SVD(X,D,U,V);                              // X = U * D * V.t()
  201.    ColumnVector Fitted = U.t() * Y;
  202.    ColumnVector A = V * ( D.i() * Fitted );
  203.    Fitted = U * Fitted;
  204.    ColumnVector Residual = Y - Fitted;
  205.    Real ResVar = Residual.SumSquare() / (nobs-npred1);
  206.  
  207.    // get variances of estimates
  208.    D << V * (D * D).i() * V.t();
  209.  
  210.    // Get diagonals of Hat matrix
  211.    DiagonalMatrix Hat;  Hat << U * U.t();
  212.  
  213.    // print out answers
  214.    cout << "\nEstimates and their standard errors\n\n";
  215.    ColumnVector SE(npred1);
  216.    for (int i=1; i<=npred1; i++) SE(i) = sqrt(D(i)*ResVar);
  217.    cout << setw(11) << setprecision(5) << (A | SE) << endl;
  218.    cout << "\nObservations, fitted value, residual value, hat value\n";
  219.    cout << setw(9) << setprecision(3) << 
  220.       (X.Columns(2,3) | Y | Fitted | Residual | Hat.AsColumn());
  221.    cout << "\n\n";
  222. }
  223.  
  224. main()
  225. {
  226.    cout << "\nDemonstration of Matrix package\n";
  227.  
  228.    // Test for any memory not deallocated after running this program
  229.    Real* s1; { ColumnVector A(8000); s1 = A.Store(); }
  230.  
  231.    {
  232.       // the data
  233.  
  234. #ifndef ATandT
  235.       Real y[]  = { 8.3, 5.5, 8.0, 8.5, 5.7, 4.4, 6.3, 7.9, 9.1 };
  236.       Real x1[] = { 2.4, 1.8, 2.4, 3.0, 2.0, 1.2, 2.0, 2.7, 3.6 };
  237.       Real x2[] = { 1.7, 0.9, 1.6, 1.9, 0.5, 0.6, 1.1, 1.0, 0.5 };
  238. #else             // for compilers that don't understand aggregrates
  239.       Real y[9], x1[9], x2[9];
  240.       y[0]=8.3; y[1]=5.5; y[2]=8.0; y[3]=8.5; y[4]=5.7;
  241.       y[5]=4.4; y[6]=6.3; y[7]=7.9; y[8]=9.1;
  242.       x1[0]=2.4; x1[1]=1.8; x1[2]=2.4; x1[3]=3.0; x1[4]=2.0;
  243.       x1[5]=1.2; x1[6]=2.0; x1[7]=2.7; x1[8]=3.6;
  244.       x2[0]=1.7; x2[1]=0.9; x2[2]=1.6; x2[3]=1.9; x2[4]=0.5;
  245.       x2[5]=0.6; x2[6]=1.1; x2[7]=1.0; x2[8]=0.5;
  246. #endif
  247.  
  248.       int nobs = 9;                           // number of observations
  249.       int npred = 2;                          // number of predictor values
  250.  
  251.       // we want to find the values of a,a1,a2 to give the best
  252.       // fit of y[i] with a0 + a1*x1[i] + a2*x2[i]
  253.       // Also print diagonal elements of hat matrix, X*(X.t()*X).i()*X.t()
  254.  
  255.       // this example demonstrates four methods of calculation
  256.  
  257.       Try
  258.       {
  259.          test1(y, x1, x2, nobs, npred);
  260.          test2(y, x1, x2, nobs, npred);
  261.          test3(y, x1, x2, nobs, npred);
  262.          test4(y, x1, x2, nobs, npred);
  263.       }
  264.       Catch(DataException) { cout << "\nInvalid data\n"; }
  265.       Catch(SpaceException) { cout << "\nMemory exhausted\n"; }
  266.       CatchAll { cout << "\nUnexpected program failure\n"; }
  267.    }
  268.  
  269. #ifdef DO_FREE_CHECK
  270.    FreeCheck::Status();
  271. #endif
  272.    Real* s2; { ColumnVector A(8000); s2 = A.Store(); }
  273.    cout << "\n\nChecking for lost memory: "
  274.       << (unsigned long)s1 << " " << (unsigned long)s2 << " ";
  275.    if (s1 != s2) cout << " - error\n"; else cout << " - ok\n";
  276.  
  277.    return 0;
  278.  
  279. }
  280.  
  281.